rm(list=ls())
library(tidyverse)
library(weights)
library(questionr)
`%notin%` = Negate(`%in%`)

setwd("~/Downloads/dataverse_files(3)")

public <- read_csv("Data/dfp_furnas_micro - dfp_furnas_micro.csv")


respondents <- read_csv("Data/CleanedDeIdentifiedSPEPS02022021_Science.csv")


respondents <- respondents %>% mutate(pid = case_when(party.main == "a Republican" &  party.rep.probe == "Strong" ~ "Strong Rep",
                                                      party.main == "a Republican" &  party.rep.probe == "Not very strong" ~ "Weak Rep",
                                                      party.main == "a Democrat" &  party.dem.probe == "Strong" ~ "Strong Dem",
                                                      party.main == "a Democrat" &  party.dem.probe == "Not very strong" ~ "Weak Dem",
                                                      party.main != "a Republican" &  party.main != "a Democrat" & party.lean.probe == "Democratic" ~ "Dem Lean",
                                                      party.main != "a Republican" &  party.main != "a Democrat" & party.lean.probe == "Republican" ~ "Rep Lean",
                                                      party.main != "a Republican" &  party.main != "a Democrat" & party.lean.probe == "Neither" ~ "Independent"
))
respondents <- respondents %>% mutate(pid3 = case_when(grepl("Dem",pid, fixed=T) ~ "Dem", 
                                                       grepl("Rep",pid, fixed=T) ~ "Rep",
                                                       grepl("Independent",pid, fixed=T) ~ "Ind",
                                                       TRUE ~ "DK"))

respondents$pid3[respondents$pid3 == "DK"] <- "Ind"
respondents$pid3[respondents$pid3 == "Ind"] <- "Ind/DK"


sum(respondents$weight_use)


round(wtd.chi.sq(respondents$trust_1_13, respondents$pid3, weight = respondents$weight_use), digits =5)


round(wtd.chi.sq(respondents$trust_1_14, respondents$pid3, weight = respondents$weight_use), digits =5)



respondents$trust_1_13 <- fct_relevel(respondents$trust_1_13,"Very untrustworthy", "Somewhat untrustworthy", "Neither", "Somewhat trustworthy", "Very trustworthy")
round(prop.table(wtd.table(respondents$trust_1_13, respondents$pid3, weights = respondents$weight_use), margin = 2)*100, digits = 1)

respondents$trust_1_14 <- fct_relevel(respondents$trust_1_14,"Very untrustworthy", "Somewhat untrustworthy", "Neither", "Somewhat trustworthy", "Very trustworthy")
round(prop.table(wtd.table(respondents$trust_1_14, respondents$pid3, weights = respondents$weight_use), margin = 2)*100, digits = 1)

respondents_sci_trust <- respondents %>% dplyr::select( trust_1_13, trust_1_14, pid3, weight_use)

respondents_sci_trust <- respondents_sci_trust %>% pivot_longer(cols = starts_with("trust"), names_to = "source", values_to = "trust_level")

respondents_sci_trust <- respondents_sci_trust %>% mutate(source = case_when(source == "trust_1_13" ~ "American Association for the Advancement of Science",
                                                                             source == "trust_1_14" ~ "National Academies of Science, Engineering, and Medicine" 
                                                                             ))
respondents_sci_trust$trust_level <- fct_relevel(respondents_sci_trust$trust_level,"Very untrustworthy", "Somewhat untrustworthy", "Neither", "Somewhat trustworthy", "Very trustworthy")

respondents_sci_trust_sum <- respondents_sci_trust %>% filter(!is.na(trust_level)) %>% group_by(pid3, source, trust_level) %>% summarise(cell_mass = sum(weight_use))
respondents_sci_trust_margin <- respondents_sci_trust_sum %>% group_by(pid3, source) %>% summarise(margin_sum = sum(cell_mass))
respondents_sci_trust_sum <- respondents_sci_trust_sum %>% left_join(respondents_sci_trust_margin) %>% mutate(proportion = cell_mass/margin_sum)

respondents_sci_trust_sum$pid3 <- fct_rev(respondents_sci_trust_sum$pid3)
respondents_sci_trust_sum$source <- fct_relevel(respondents_sci_trust_sum$source, "National Academies of Science, Engineering, and Medicine","American Association for the Advancement of Science")

respondents_sci_trust_sum <- respondents_sci_trust_sum %>% mutate(trust_binary = case_when(trust_level == "Very untrustworthy" ~ "Untrustworthy", 
                                                                                           trust_level == "Somewhat untrustworthy" ~ "Untrustworthy",
                                                                                           trust_level == "Very trustworthy" ~ "Trustworthy", 
                                                                                           trust_level == "Somewhat trustworthy" ~ "Trustworthy",
                                                                                           trust_level == "Neither" ~ "Neither"))

respondents_sci_trust_sum <- respondents_sci_trust_sum %>% mutate(very_binary = case_when(trust_level == "Very untrustworthy" ~ "Very", 
                                                                                           trust_level == "Somewhat untrustworthy" ~ "Somewhat",
                                                                                           trust_level == "Very trustworthy" ~ "Very", 
                                                                                           trust_level == "Somewhat trustworthy" ~ "Somewhat",
                                                                                           trust_level == "Neither" ~ "Neither"))
respondents_sci_trust_sum$pid3 <- respondents_sci_trust_sum$pid3 %>% fct_rev()
respondents_sci_trust_sum$very_binary <- respondents_sci_trust_sum$very_binary %>% fct_rev()

respondents_sci_trust_sum$color_group <- paste(respondents_sci_trust_sum$pid3, respondents_sci_trust_sum$trust_level, sep = " - ")

respondents_sci_trust_sum <- respondents_sci_trust_sum %>% mutate(source_label = case_when(source == "National Academies of Science, Engineering, and Medicine" ~ "National Academies of Science,\nEngineering, and Medicine", 
                                                                                           source ==  "American Association for the Advancement of Science" ~ "American Association for the\nAdvancement of Science",
                                                                                          TRUE ~ as.character(source)))
respondents_sci_trust_sum$source_label <- fct_relevel(respondents_sci_trust_sum$source_label, "National Academies of Science,\nEngineering, and Medicine","American Association for the\nAdvancement of Science","Union of Concerned Scientists"  )

(p1_use <- ggplot(filter(respondents_sci_trust_sum, trust_binary != "Neither", trust_binary != "Untrustworthy") , aes(x=pid3, y= proportion, group= trust_level, fill = color_group)) + 
    geom_col(width=.5) +
    scale_fill_manual("", values=c("#2166ac","#053061","#878787","#666666", "#b2182b", "#67001f"), guide = guide_legend(reverse = F)) + facet_grid(trust_binary~source_label) +
    theme_minimal() + scale_y_continuous(labels = scales::percent) + ylab("Proportion of Political Elites")  +
    theme(legend.position = "left", strip.text.x = element_text(size = 12,face = "bold", hjust = 0), axis.text.x = element_text(color = c("#2166ac","#666666", "#b2182b"), face= "bold", size=12)) +
    xlab("")
)   
    
ggsave("Output/FigS27.pdf",  p1_use, device = "pdf", width = 5, height = 3, units="in")









library(cowplot)

round(wtd.chi.sq(respondents$trust_science_1, respondents$pid3, weight = respondents$weight_use), digits =5)
round(wtd.chi.sq(respondents$trust_science_2, respondents$pid3, weight = respondents$weight_use), digits =5)
round(wtd.chi.sq(respondents$trust_science_3, respondents$pid3, weight = respondents$weight_use), digits =5)
round(wtd.chi.sq(respondents$trust_science_4, respondents$pid3, weight = respondents$weight_use), digits =5)


respondents$trust_science_num_1 <- fct_relevel(respondents$trust_science_1, "Completely distrust", "Partially distrust", "Neither distrust nor trust", "Partially trust", "Completely trust") %>% as.numeric() - 3
respondents$trust_science_num_2 <- fct_relevel(respondents$trust_science_2, "Completely distrust", "Partially distrust", "Neither distrust nor trust", "Partially trust", "Completely trust") %>% as.numeric() - 3
respondents$trust_science_num_3 <- fct_relevel(respondents$trust_science_3, "Completely distrust", "Partially distrust", "Neither distrust nor trust", "Partially trust", "Completely trust") %>% as.numeric() - 3
respondents$trust_science_num_4 <- fct_relevel(respondents$trust_science_4, "Completely distrust", "Partially distrust", "Neither distrust nor trust", "Partially trust", "Completely trust") %>% as.numeric() - 3

public$dfp_furnas_trust_science_num_1 <- fct_relevel(public$dfp_furnas_trust_science_1 , "Completely distrust", "Partially distrust", "Neither distrust or trust", "Partially trust", "Completely trust") %>% as.numeric() - 3
public$dfp_furnas_trust_science_num_2 <- fct_relevel(public$dfp_furnas_trust_science_2, "Completely distrust", "Partially distrust", "Neither distrust or trust", "Partially trust", "Completely trust") %>% as.numeric() - 3
public$dfp_furnas_trust_science_num_3 <- fct_relevel(public$dfp_furnas_trust_science_3, "Completely distrust", "Partially distrust", "Neither distrust or trust", "Partially trust", "Completely trust") %>% as.numeric() - 3
public$dfp_furnas_trust_science_num_4 <- fct_relevel(public$dfp_furnas_trust_science_4, "Completely distrust", "Partially distrust", "Neither distrust or trust", "Partially trust", "Completely trust") %>% as.numeric() - 3


respondents <- respondents %>% mutate(trust_science_mean = rowMeans(select(respondents, starts_with("trust_science_num")), na.rm = TRUE))
respondents$trust_science_mean

library(survey)

respondents_svy<-svydesign(id = ~1, weights = ~weight_use, data = respondents)


wtd.t.test(x=respondents$trust_science_mean[respondents$pid3 == "Dem"], y = respondents$trust_science_mean[respondents$pid3 == "Rep"], weight =respondents$weight_use[respondents$pid3 == "Dem"], weighty =respondents$weight_use[respondents$pid3 == "Rep"], bootse = T, bootp=T, bootn=10000)



respondents_scientist_trust <- respondents %>% dplyr::select(pid3, weight_use, trust_science_1, trust_science_2, trust_science_3, trust_science_4)

public_scientist_trust <- public %>% dplyr::select(partisanship, weight_column, dfp_furnas_trust_science_1, dfp_furnas_trust_science_2, dfp_furnas_trust_science_3, dfp_furnas_trust_science_4)
public_scientist_trust <- public_scientist_trust %>% mutate(partisanship =case_when(partisanship == "Democrat" ~ "Dem", 
                                                                                    partisanship == "Republican" ~ "Rep",
                                                                                    partisanship == "Independent / Third Party" ~ "Ind/DK"))


public_scientist_trust[public_scientist_trust == "Neither distrust or trust"] <- "Neither distrust nor trust"   



names(public_scientist_trust) <- c("pid3", "weight_use", "trust_science_1", "trust_science_2", "trust_science_3", "trust_science_4")
public_scientist_trust$sample <- "Likely Voters"
respondents_scientist_trust$sample <- "Political Elites"
respondents_scientist_trust <- bind_rows(respondents_scientist_trust, public_scientist_trust)

respondents_scientist_trust_long <- respondents_scientist_trust %>% pivot_longer(cols = starts_with("trust"), names_to = "prompt", values_to = "trust_level")

respondents_scientist_trust_long$trust_level <- fct_relevel(respondents_scientist_trust_long$trust_level, "Completely distrust", "Partially distrust", "Neither distrust nor trust", "Partially trust", "Completely trust")

respondents_scientist_trust_sum <- respondents_scientist_trust_long %>% filter(!is.na(trust_level)) %>% group_by(pid3, prompt, trust_level, sample) %>% summarise(cell_mass = sum(weight_use))
respondents_scientist_trust_margin <- respondents_scientist_trust_sum %>% group_by(pid3, prompt, sample) %>% summarise(margin_sum = sum(cell_mass))
respondents_scientist_trust_sum <- respondents_scientist_trust_sum %>% left_join(respondents_scientist_trust_margin) %>% mutate(proportion = cell_mass/margin_sum)

respondents_scientist_trust_sum$pid3 <- fct_rev(respondents_scientist_trust_sum$pid3)
respondents_scientist_trust_sum$trust_level <- fct_relevel(respondents_scientist_trust_sum$trust_level, "Completely distrust", "Partially distrust", "Neither distrust nor trust", "Partially trust", "Completely trust")

respondents_scientist_trust_sum <- respondents_scientist_trust_sum %>% mutate(prompt_long = case_when(prompt == "trust_science_1" ~ "Create knowledge that is unbiased and accurate", 
                                                                                                      prompt == "trust_science_2" ~ "Create knowledge that is useful",
                                                                                                      prompt == "trust_science_3" ~ "Advise government officials on policy",
                                                                                                      prompt == "trust_science_4" ~ "Inform the public on important issues"
))

respondents_scientist_trust_sum <- respondents_scientist_trust_sum %>% mutate(trust_binary = case_when(trust_level == "Completely distrust" ~ "Distrust", 
                                                                                           trust_level == "Partially distrust" ~ "Distrust",
                                                                                           trust_level == "Completely trust" ~ "Trust", 
                                                                                           trust_level == "Partially trust" ~ "Trust",
                                                                                           trust_level == "Neither" ~ "Neither"))

respondents_scientist_trust_sum <- respondents_scientist_trust_sum %>% mutate(complete_binary = case_when(trust_level == "Completely distrust" ~ "Complete", 
                                                                                                      trust_level == "Partially distrust" ~ "Partially",
                                                                                                      trust_level == "Completely trust" ~ "Complete", 
                                                                                                      trust_level == "Partially trust" ~ "Partially",
                                                                                                      trust_level == "Neither" ~ "Neither"))





respondents_scientist_trust_sum$pid3 <- respondents_scientist_trust_sum$pid3 %>% fct_rev()
respondents_scientist_trust_sum$complete_binary <- respondents_scientist_trust_sum$complete_binary %>% fct_rev()
respondents_scientist_trust_sum$trust_binary <- respondents_scientist_trust_sum$trust_binary %>% fct_rev()

respondents_scientist_trust_sum$color_group <- paste(respondents_scientist_trust_sum$pid3, respondents_scientist_trust_sum$trust_level, sep = " - ")


(p3.3 <- ggplot(filter(respondents_scientist_trust_sum,trust_level %in% c("Completely trust","Partially trust" )) , aes(x=pid3, y= proportion, group= trust_level, fill = color_group)) + 
    geom_col(width=.4) + 
    scale_fill_manual("", values=c("#053061","#2166ac","#666666","#878787", "#67001f","#b2182b"), guide = guide_legend(reverse = TRUE)) + facet_grid(sample~prompt_long) +
    theme_minimal() + scale_y_continuous(labels = scales::percent) + ylab("Proportion of Political Elites")  + xlab("")+
    theme(legend.position = "bottom", strip.text.x = element_text(size = 12,face = "bold", hjust = 0), axis.text.x = element_text(color = c("#2166ac","#666666", "#b2182b"), face= "bold", size=12))
) 
ggsave("Output/Fig4.pdf",  p3.3, device = "pdf", width = 8, height = 6, units="in")


respondents$racial_resentment_num_1 <- fct_relevel(respondents$racial_resentment_1, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree")%>% as.numeric() -3
respondents$racial_resentment_num_2 <- fct_relevel(respondents$racial_resentment_2, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree")%>% as.numeric() -3
respondents$racial_resentment_num_3 <- fct_relevel(respondents$racial_resentment_3, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree")%>% fct_rev() %>% as.numeric() -3
respondents <- respondents %>% mutate(racial_resentment_mean = rowMeans(select(respondents, starts_with("racial_resentment_num")), na.rm = TRUE))


respondents$sexism_num_1 <- fct_relevel(respondents$sexism_1, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree") %>% fct_rev()  %>% as.numeric() -3
respondents$sexism_num_2 <- fct_relevel(respondents$sexism_2, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree") %>% fct_rev() %>% as.numeric() -3
respondents$sexism_num_3 <- fct_relevel(respondents$sexism_3, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree") %>% fct_rev() %>% as.numeric() -3
respondents <- respondents %>% mutate(sexism_mean = rowMeans(select(respondents, starts_with("sexism_num")), na.rm = TRUE))

respondents <- respondents %>% mutate(conspiracy_mean = rowMeans(select(respondents, conspiracy_1, conspiracy_2, conspiracy_3), na.rm = TRUE))

respondents$trust_1_13 <- fct_relevel(respondents$trust_1_13,"Very untrustworthy", "Somewhat untrustworthy", "Neither", "Somewhat trustworthy", "Very trustworthy")  %>% as.numeric() -3

respondents$trust_1_14 <- fct_relevel(respondents$trust_1_14,"Very untrustworthy", "Somewhat untrustworthy", "Neither", "Somewhat trustworthy", "Very trustworthy")  %>% as.numeric() -3

respondents <- respondents %>% mutate(science_institutions_trust_mean = rowMeans(select(respondents, c(trust_1_13, trust_1_14)), na.rm = TRUE))


library(eRm)
library(forcats)
respondents$ideo.bat_1_num <- fct_relevel(respondents$ideo.bat_1, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree") %>% fct_rev() %>% as.numeric()
respondents$ideo.bat_2_num <- fct_relevel(respondents$ideo.bat_2, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree") %>% as.numeric()
respondents$ideo.bat_3_num <- fct_relevel(respondents$ideo.bat_3, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree")%>% as.numeric()
respondents$ideo.bat_4_num <- fct_relevel(respondents$ideo.bat_4, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree")%>% as.numeric()
respondents$ideo.bat_5_num <- fct_relevel(respondents$ideo.bat_5, "Strongly agree", "Somewhat agree", "Neither agree nor disagree", "Somewhat disagree", "Strongly disagree")%>% as.numeric()



dat_ideo <- respondents %>% dplyr::select(ResponseId, ideo.bat_1_num, ideo.bat_2_num,ideo.bat_3_num,ideo.bat_4_num,ideo.bat_5_num)
dat_ideo <- dat_ideo %>% na.omit


res.rsm <- PCM(dplyr::select(dat_ideo, -ResponseId))
thresholds(res.rsm)
par(mfrow=c(3,2))
plotICC(res.rsm, mplot = FALSE, legpos = FALSE, ask = FALSE)

plotPImap(res.rsm, sorted = TRUE)
pf <- person.parameter(res.rsm)
ifs <- itemfit(pf)
ifs

library(weights)
dat_ideo$ideology.irt <- coef(pf)
respondents <- respondents %>% left_join(dplyr::select(dat_ideo, ResponseId, ideology.irt))
respondents <- respondents %>% mutate(ideology.irt.std = stdz(ideology.irt))



respondents$pid <- fct_relevel(respondents$pid, "Independent", "Strong Dem", "Weak Dem", "Dem Lean", "Rep Lean", "Weak Rep", "Strong Rep")
m_ts <- lm(stdz(trust_science_mean, weight_use) ~ pid  + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = respondents)
m_tsi <- lm(stdz(science_institutions_trust_mean, weight_use) ~ pid  + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = respondents)

library(modelsummary)

modelsummary(list("Trust in scientists" = m_ts, "Trust in Scientific Institutions" = m_tsi), output = "Output/TableS31.docx")

